#include "libmri.h"

static void iSPIRiT (cplxf *y, bool adj, 
	cplxf *k, int ks1, int ks2, int ch1, int ch2,
	cplxf *x, int im1, int im2, int ch);

static void kSPIRiT (cplxf *y, bool adj, 
	cplxf *k, int ks1, int ks2, int ch1, int ch2,
	cplxf *x, int im1, int im2, int ch);

extern void spirit (cplxf *y,
	bool adj,
	cplxf *k, int ks1, int ks2, int kch1, int kch2,
	cplxf *x, int im1, int im2, int xch)
{
	if (adj){ // assert xch == kch2
		for (int k1 = 0; k1 < kch1; ++k1)
			for (int j = 0; j < im2; ++j)
				for (int i = 0; i < im1; ++i)
					y[i + im1*(j + im2*k1)] = 0.0f;
	}
	else { // assert xch == kch1
		for (int k2 = 0; k2 < kch2; ++k2)
			for (int j = 0; j < im2; ++j)
				for (int i = 0; i < im1; ++i)
					y[i + im1*(j + im2*k2)] = 0.0f;
	}

	if (ks1 == im1 && ks2 == im2){
		iSPIRiT (y, adj, k, ks1, ks2, kch1, kch2, x, im1, im2, xch);
	}
	else{
		kSPIRiT (y, adj, k, ks1, ks2, kch1, kch2, x, im1, im2, xch);
	}
}

static void iSPIRiT (cplxf *y, bool adj, 
	cplxf *k, int ks1, int ks2, int kch1, int kch2,
	cplxf *x, int im1, int im2, int xch)
{
	if (adj){ // assert xch == kch2
		for (int k1 = 0; k1 < kch1; ++k1)
			for (int k2 = 0; k2 < kch2; ++k2)
				for (int j = 0; j < im2; ++j)
					for (int i = 0; i < im1; ++i)
					{
						y[i + im1*(j + im2*k1)]
							+= x[i + im1*(j + im2*k2)]
							   * conj (k[i + im1*(j + im2*(k1 + kch1*k2))]);
					}
	}
	else { // assert xch == kch1
		for (int k2 = 0; k2 < kch2; ++k2)
			for (int k1 = 0; k1 < kch1; ++k1)
				for (int j = 0; j < im2; ++j)
					for (int i = 0; i < im1; ++i)
					{
						y[i + im1*(j + im2*k2)]
							+= x[i + im1*(j + im2*k1)]
							   * k[i + im1*(j + im2*(k1 + kch1*k2))];
					}
	}
}

static void kSPIRiT (cplxf *y, bool adj, 
	cplxf *k, int ks1, int ks2, int kch1, int kch2,
	cplxf *x, int im1, int im2, int xch)
{
	if (adj){ // assert xch == kch2
		for (int c1 = 0; c1 < kch1; ++c1)
			for (int c2 = 0; c2 < kch2; ++c2)
			{
				cplxf *gij = &k[ks1*ks2*(c1 + kch1*c2)];

				for (int j = 0; j < im2; ++j)
					for (int i = 0; i < im1; ++i)
					{
						for (int k2 = 0; k2 < ks2; ++k2)
							for (int k1 = 0; k1 < ks1; ++k1)
							{
								int xi = i-(ks1/2)+k1;
								int xj = j-(ks2/2)+k2;

								if (xi >= 0 && xi < im1 && xj >= 0 && xj < im2)
									y[i + im1*(j + im2*c1)]
										+= conj(gij[ks1-1-k1 + ks1*(ks2-1-k2)])
										   * x[xi + im1*(xj + im2*c2)];
							}
					}
			}
	}
	else { // assert xch == kch1
		for (int c2 = 0; c2 < kch2; ++c2)
			for (int c1 = 0; c1 < kch1; ++c1)
			{
				for (int j = 0; j < im2; ++j)
					for (int i = 0; i < im1; ++i)
					{
						for (int k2 = 0; k2 < ks2; ++k2)
							for (int k1 = 0; k1 < ks1; ++k1)
							{
								int xi = i-(ks1/2)+k1;
								int xj = j-(ks2/2)+k2;

								if (xi >= 0 && xi < im1 && xj >= 0 && xj < im2)
									y[i + im1*(j + im2*c2)] += 
										k[k1 + ks1*(k2 + ks2*(c1 + kch1*c2))]
										* x[xi + im1*(xj + im2*c1)];
							}
					}
			}
	}
}
